clear
rng default

%% Load useful matrices
load PYX_cond 
load U_CS
x=3;
PYX_cond_x=PYX_cond(x,:);

%% Useful parameters 
n=8^5;
p=n*(n-1)/2;
m_R1=8*10^6;
random_indices_pairs_original_R1 = randperm(p, m_R1).';

n=15^5;
p=n*(n-1)/2;
m_R2=4*10^6;
random_indices_pairs_original_R2 = randperm(p, m_R2).';

n=8^5;
p=n*(n-1)/2;
m_restricted_all_R1=[10^5 m_R1]; 
random_indices_pairs_all_R1=cell(size(m_restricted_all_R1,2),1);
random_indices_pairs_all_R1{end}=random_indices_pairs_original_R1;
for t=1:size(m_restricted_all_R1,2)-1
    random_indices_pairs_all_R1{t}=randperm(p, m_restricted_all_R1(t)).';
end

n=15^5;
p=n*(n-1)/2;
m_restricted_all_R2=[10^5 m_R2]; 
random_indices_pairs_all_R2=cell(size(m_restricted_all_R2,2),1);
random_indices_pairs_all_R2{end}=random_indices_pairs_original_R2;
for t=1:size(m_restricted_all_R2,2)-1
    random_indices_pairs_all_R2{t}=randperm(p, m_restricted_all_R2(t)).';
end

clear random_indices_pairs_original_R1 random_indices_pairs_original_R2

Ueq=0;
Uineq=0;
tol=10^(-4);

param_MOSEK.MSK_IPAR_LOG = 0; 
param_MOSEK.MSK_IPAR_INTPNT_BASIS = 'MSK_BI_NEVER';
param_MOSEK.MSK_IPAR_NUM_THREADS = 1;

%% Construct grid
u0grid=0;
u1grid=U_CS(3,1); 
gran=60; 
% try different values of gran and different sizes of grid hypercube as discussed in Appendix C of the paper (need powerful cluster to run the code for higher values of gran and larger sizes of hypercube)
u2grid=linspace(-15,15,gran); 
u3grid=linspace(-15,15,gran); 
u4grid=linspace(-15,15,gran); 
u5grid=linspace(-15,15,gran); 
[ca, cb, cc,cd, ce, cf] = ndgrid(u0grid, u1grid, u2grid, u3grid, u4grid, u5grid);
u0grid=ca(:);
u1grid=cb(:);
u2grid=cc(:);
u3grid=cd(:);
u4grid=ce(:);
u5grid=ce(:);

%% Workers
workers=600; %number of parallel workers
jobs=round(size(u2grid,1)/workers); %number of jobs per parallel worker

